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PART  1:  STATISTICAL  APPROACHES  TO  SEISMIC  DISCRIMINATION 


1.1  General  Approaches  to  Discrimination 

Conventional  methods  for  discriminating  between  earthquakes  and  explosions  at  re¬ 
gional  distances  have  concentrated  on  extracting  specific  features  from  the  waveforms  of 
the  two  arrival  phases  seen  on  the  seismic  record,  denoted  in  this  report  by  P  and  S.  The 
specific  features  usually  considered  are  amplitude  ratios,  measures  of  waveform  complex¬ 
ity  or  various  kinds  of  spectral  ratios,  suggesting  that  the  main  characterization  of  the 
differences  between  earthquakes  and  explosions  can  be  reduced  to  comparing  differences 
between  the  amplitude  ratios  or  the  spectral  ratios  of  the  two  arrival  phases.  For  purposes 
of  discussion,  we  will  term  this  approach  the  feature  extraction  approach.  The  variations 
on  the  above  theme  seem  to  change  depending  on  the  author  and  on  the  general  region 
within  which  the  discrimination  study  is  done.  The  features  may  then  combined  using  an 
optimal  statistical  procedure  that  assumes  multivariate  normality  for  the  extracted  featme 
vector  and  maximizes  the  probability  of  detecting  an  explosion  for  a  fixed  specified  rate  of 
false  alarms. 

An  alternate  approach,  based  on  classical  signal  discrimination  methodology,  is  to  con¬ 
sider  the  two  P  and  S  waveforms  as  the  imderlying  phenomena  rather  than  the  somewhat 
arbitrarily  extracted  features  discussed  above.  The  two  waveforms  are  then  modeled  as 
general  zero-mean  stochastic  processes  with  given  unequal  covariance  functions  and  den¬ 
sity  functions  based  on  the  multivariate  normal  assumption.  The  same  optimality  criterion 
as  above  is  applied  to  the  underlying  waveforms.  We  shall  term  this  the  waveform  dis¬ 
crimination  approach.  This  approach  has  the  advantage  (or  disadvantage,  depending  on 
one’s  point  of  view)  of  eliminating  the  arbitrary  step  that  decides  in  advance  that  certain 
extracted  features  are  optimal  and  lets  the  statistical  model  dictate  the  method  for  com¬ 
bining  the  elements  of  the  underlying  signal  waveforms.  Modifications  that  are  robust  to 
departures  from  the  above  assiunptions  can  be  made  using  information  theoretic  measures. 

It  is  natural  and  expected  that  seismologists  are  more  comfortable  with  the  feature 
extraction  approach  (Blandford,  1993,  Dysart  and  Pulli,  1990,  Fisk  et  al,  1993,  Richards 
et  al,  1993,  Ryall,  1993,  Taylor  et  al,  1989,  Booker  and  Mitronovas,  1964)  and  that  en¬ 
gineers  and  statisticians,  who  may  lack  the  requisite  geophysical  knowledge  and  intuition 
for  extracting  reasonable  features,  should  prefer  the  signal  discrimination  approach  (Ca¬ 
vanaugh  et  al,  1993,  Dargahi-Noubary  and  Laycock,  1981,  Liggett,  1971,  Shumway,  1982, 
1983,  Shmnway  and  Unger,  1974,  Shumway  and  Blandford,  1974,  Zhang  and  Taniguchi, 
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1993,  1994).  Happily,  the  final  results  of  both  procedures  depend  upon  the  spectra  of  the 
two  phases  and  the  two  methods  seem  to  yield  somewhat  similar  results  when  applied  to 
sample  earthquakes  and  explosions. 

The  organization  of  Part  1  of  this  report  is  to  first  discuss  the  statistical  approaches  to 
discrimination  that  hold  regardless  of  whether  the  underlying  elements  are  the  extracted 
featmes  or  the  root  waveforms.  This  is  followed  by  a  more  specific  discussion  of  the  two 
approaches  with  applications  to  a  population  of  eight  earthquakes,  eight  mining  explosions 
in  Scandinavia  and  the  Novaya  Zemlya  event  (see  Ryall,  1993)  of  unknown  origin.  The  per¬ 
formance  of  the  two  approaches  is  compared  on  the  common  events  and  for  the  new  event 
of  unknown  origin;  it  is  shown  that  a  robust  modification  to  the  waveform  discriminator 
performs  best  in  this  particular  situation. 

As  demonstration  data,  we  use  a  subset  of  stations  recording  8  earthquakes  and  8 
explosions  in  Scandinavia  from  the  arrays  NORESS,  ARCESS  and  FINESS  as  described 
in  Blandford  [4]  (see  also  Cavanaugh  et  al  (1993)  for  more  details).  According  to  Blandford, 
“The  events  were  selected  with  consideration  for  having  sufficient  S/N  at  single  elements 
so  that  all  phases  could  be  clearly  seen  on  all  components  of  a  single  instrument  •••”. 
All  events  chosen  by  Blandford  were  on  or  near  land  and  were  distributed  uniformly  over 
Scandinavia  to  minimize  the  possibility  that  discriminators  might  be  keying  on  location 
or  land-sea  differences.  Figure  1  shows  portions  of  a  typical  earthquake  and  explosion 
(sampled  at  20  Hz),  along  with  the  unknown  Novaya  Zemyla  event  of  31  December,  1992. 
The  earthquake  exhibits  the  often  observed  lower  amplitude  of  the  first  arrival  P  phase  as 
compared  to  the  later  S  phase.  The  explosions  ranged  from  magnitude  2.13  to  2.19  with 
the  earthquakes  in  the  range  2.74  to  4.40.  The  Novaya  Zemyla  event  was  of  magnitude  2.5. 
For  computations,  we  did  not  identify  specific  phases  through  velocity  computations  but 
simply  chose  fairly  broad  (25  second)  windows  that  seemed  to  include  the  major  portions 
of  the  P  and  S  phases. 

1.2  Statistical  Optimality 

The  problem  of  discriminating  between  earthquakes  and  explosions  is  an  example  of 
a  general  pattern  recognition  problem  in  statistics  where  we  wish  to  classify  an  observed 
vector  X  =  (xi,X2, . . .  iXp)'  consisting  of  p  features  or  of  a  p  dimensional  waveform  into 
one  of  two  populations,  denoted  in  general  by  Hi:  Earthquake  and  H2:  Explosion.  Key 
statistical  properties  of  any  classification  procedure  are  the  false  alarm  rate  P(2|l),  defined 
as  the  probability  of  accepting  the  explosion  hypothesis,  given  that  the  event  is  an 
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earthqucike  and  P(2|2),  the  probability  of  detecting  an  explosion.  The  seismic  discrimina¬ 
tion  problem  is  expressed  then  as  one  of  distinguishing  between  two  hypotheses  and  fits 
nicely  into  classical  signal  discrimination  and  detection  theory.  Detailed  treatments  can 
be  found  in  many  places;  we  mention  Anderson  (1971)  for  an  exposition  of  the  background 
statistics  and  Shumway  (1982),  (1988)  for  extensions  of  these  ideas  to  time  series  analysis. 

To  begin,  suppose  that  we  denote  the  probability  densities  of  the  observed  vector 
under  the  two  hypotheses  by  Pi(x)  and  P2(x)  respectively.  It  is  well  known  that  the  best 
rule,  in  the  sense  of  maximizing  the  detection  probability  for  a  fixed  false  alarm  rate, 
results  from  classifying  x  into  Hi  when 

EiW>a' 

P2(x)  ^ 

and  into  H2  otherwise.  Simple  rules  result  when  one  further  specializes  by  requiring  that 
the  densities  pi(x)  and  P2(x)  correspond  to  multivariate  normal  distributions.  One  chooses 
the  constant  to  have  a  specified  false  alarm  rate  or  as  the  ratio  of  prior  probabilities  of 
H2  and  Hi  in  the  Bayesian  case.  It  is  plausible  to  assume  that  the  two  populations  differ 
only  in  the  mean  vectors  when  one  is  looking  at  extracted  features  and  differ  only  in  the 
covariance  structure  when  one  is  looking  at  waveforms. 

For  the  case  where  only  the  means  fj,i  and  ^2  are  different  and  the  covariance  matrices 
Ri  and  R2  take  some  common  value  R,  the  rule  given  above  reduces  to  choosing  Hi 
whenever  C2  >  C2  and  choosing  H2  otherwise,  where 

Cj  =  fx'j  R-^x  -  ^  n'j  (1-1) 

j  =  1,2,  and  we  assume  K  =  1  corresponding  to  equal  prior  probabilities  for  Hi  and 
H2  ■  This  linear  discriminant  function  is  easy  to  apply  and  leads  to  comparing  two  single 
numbers  in  order  to  make  a  decision.  If  the  mean  vectors  and  covariance  matrix  are 
known,  the  distributional  properties  of  the  linear  function  can  be  evaluated  explicitly  to 
give  predicted  values  for  the  performance  measures  P(2|l)  and  P(2|2).  We  temporarily 
defer  discussion  of  what  to  do  when  this  is  not  the  case. 

When  the  means  /ii  =  /i2  =  0  are  the  same  but  the  covariance  functions  Ri  and  R2 
are  different,  the  test  procedure  reduces  to  accepting  Hi  when  di  >  6,2  and  H2  otherwise, 
where 

=  (1.2) 

j  =  1,2,  which  shows  that  the  criterion  in  this  case  reduces  to  a  quadratic  discriminant 
function.  Even  if  the  covariance  matrices  are  known,  the  distribution  of  the  test  statistic 
is  intractable  and  the  performance  measures  are  extremely  difficult  to  evaluate. 
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Both  procedures  given  above  require  estimators  for  the  means  and  covariance  matrices 
if  they  cannot  assumed  to  be  known.  In  the  seismic  case,  this  means  that  small  training 
samples,  often  highly  dependent  on  source  and  array  parameters,  must  be  used  in  order  to 
develop  sample  mean  vectors  and  covariance  matrices.  The  performance  characteristics  of 
the  linear  and  quadratic  discriminant  functions  become  even  more  uncertain  under  these 
conditions  and  we  need  to  use  resampling  techniques  to  estimate  the  performance  mea¬ 
sures  P(2|l)  and  P(2|2).  A  commonly  used  resampling  technique  that  produces  reasonable 
estimators  is  to  compute  the  holdout  classification  function  proposed  by  Lachenbruch  and 
Mickey  (1968).  In  this  procedure,  one  “holds  out”  the  observation  to  be  classified  when 
estimating  the  discriminant  functions  (1.1)  and  (1.2).  That  is,  the  sample  mean  and 
covariance  matrices  are  computed  for  the  training  sample  with  the  observation  to  be  clas¬ 
sified  held  out.  The  overall  sample  error  rates  and  detection  probabilities  are  then  better 
estimators  for  the  true  population  values. 

An  alternative  but  closely  related  proposal  by  Fisk  et  al  (1993)  derives  the  likelihood 
ratio  criterion  for  testing  whether  the  held  out  observation  belongs  to  the  earthquake 
or  explosion  population  (see  also  Anderson,  1971).  The  performance  measures  are  then 
evaluated  by  resampling  the  likelihood  ratio  test  statistic  using  the  bootstrap.  The  authors 
have  shown  in  simulations  that  this  procedure  also  produces  reasonable  estimators  for  the 
performance  measures.  We  do  not  discuss  this  interesting  proposal  further  here  and  have 
not  replicated  this  approach  using  our  data. 

It  is  clear  that  the  above  procedures  depend  more  or  less  on  a  number  of  assumptions. 
For  the  mean  differences  approach,  applied  often  in  feature  extraction,  the  computations 
are  fairly  simple  because  the  dimensionality  is  small.  The  Gaussian  assumption  is  needed 
for  the  linear  function  to  be  optimal  and  this  is  almost  always  plausible  if  one  transforms 
using  logarithms.  Equality  of  the  two  covariance  matrices  is  assumed  and  generally  can  be 
satisfied  by  using  ratios  of  features  from  the  P  and  S  phases.  The  effects  of  small  samples 
and  other  departures  are  mitigated  by  using  a  resampling  procedure. 

The  waveform  discrimination  approach  introduces  additional  computational  difficul¬ 
ties  because  of  the  large  dimensions  of  the  signal  vectors.  This  is  circumvented  by  using 
spectral  approximations  for  the  densities  pi(x)  and  P2(x)  proposed  first  by  Whittle  (1954). 
Departures  from  the  Gaussian  assumption  are  not  particularly  important  since  the  approx¬ 
imations  involve  sums  which  converge  to  Gaussian  distributions  rather  quickly  and  because 
recent  analyses  by  Zhang  and  Taniguchi  (1993),  (1994)  have  shown  that  modifications  to 
the  classical  approximations  are  robust  towards  such  departures.  Even  perturbations  in 
the  spectra  can  be  tolerated  for  the  test  statistics  used  in  Section  4.  Again,  resampling 
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techniques  are  necessary  in  order  to  provide  reasonable  estimates  for  the  false  alarm  and 
signal  detection  probabilities. 

1.3  Feature  Extraction 

In  order  to  develop  a  linear  discriminant  function  for  extracting  features  using  (1.1),  we 
need  to  determine  the  ingredients  for  the  feature  vector  x  =  (xi,  0:2, . . . , Xp)'.  An  early  ap¬ 
plication  of  this  idea  to  discriminating  at  teleseismic  distances  was  by  Booker  and  Mitrono- 
vas  (1964)  who  used  surface  wave  and  body  wave  magnitudes  as  components.  In  the  case  of 
regional  data,  numerous  investigators  have  pointed  out  that  the  logarithms  of  Pg/Lg(Pf  S) 
amplitude  ratios  tend  to  be  lower  for  earthquakes  than  for  explosions.  The  top  panel  in 
Figure  2  shows  the  values  of  —  log  P/S  amplitudes  for  the  8  earthquakes  and  explosions 
and  for  Novaya  Zemyla.  It  is  clear  that  there  is  reasonable  separation  and  that  the  Novaya 
Zemlya  event  falls  with  the  explosion  group.  Considerable  past  effort  has  been  expended 
on  spectral  ratios  involving  the  P  and  S  groups.  Bennett  and  Murphy  (1986)  note  that 
for  western  U.S.  events,  earthquake  Lg  spectra  contained  more  high  frequencies,  and  that 
the  ratio  of  the  logarithms  of  low  frequency  (.5-1  Hz)  Lg  to  higher  frequency  Lg  (2-4  Hz) 
tend  to  be  larger  for  explosions.  Taylor  et  al  (1989)  also  use  this  ratio  over  the  frequency 
bands  (1-2  Hz)  and  (6-8  Hz)  and  extended  the  consideration  to  the  Pg  phase.  Dysart  and 
PuUi  (1990)  have  also  considered  various  spectral  ratios  P/S  for  Scandinavian  events  and 
have  developed  nonlinear  neural  networks  as  an  alternative  to  simple  linear  combinations 
of  features  for  discrimination.  They  note  that  the  P/  S  spectral  ratios  are  generally  higher 
for  explosions  than  for  earthquakes.  Finally,  Richards  et  al  (1993)  note  that  for  eastern 
U.S.  events  the  ratios  of  Pg  to  Lg  spectra  are  generally  higher  for  explosions.  For  our 
particular  feature  extraction  example  using  the  Scandinavian  earthquakes  and  explosions, 
we  consider  incorporating  P  and  S  amplitudes  and  P  and  S  spectra  over  relatively  broad 
low  and  high  frequency  bands  (0-8  Hz,  6-12  Hz).  The  frequency  ranges  were  not  exactly 
comparable  to  those  used  in  the  literature  (.5-lHz,  2-4  Hz  in  Bennett  and  Murphy,  1986, 
1-2  Hz,  6-8  Hz  in  Taylor  et  al,  1989,  2-5  Hz,  5-10  Hz,  10-20  Hz  in  Dysart  and  Pulli,  1990, 
5-25(5)  Hz  in  Richards  et  al,  1993  )  but  were  chosen  by  visually  inspecting  the  separate 
spectra  and  the  average  earthquake  and  explosion  spectra.  Taking  logarithms  improves 
the  approximation  to  normality  required  for  application  of  (1.1).  We  also  used  the  ra¬ 
tios  of  P  to  5  components  (—log P/5)  in  order  to  equalize  the  covariance  matrices  of 
the  earthquake  and  explosion  populations.  Hence,  we  ended  with  three  basic  components 
representing  amplitude  ratios  and  two  spectral  ratios,  say  x  =  (xi,X2,X3)'  for  application 
in  the  linear  discriminant  function  (1).  Figure  2  shows  the  separation  for  the  populations 
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Figure  2:  Separation  Achieved  for  Earthquakes  and  Explosions  Using  Log  Amphtude  Ratios 
(Top  Panel)  Low  Frequency  (0-8Hz)  (Middle  Panel)  and  High  Frequency  (6-  12 
Hz)  (Bottom  Panel)  Spectral  Ratios. 
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of  8  earthquakes  and  explosions  spectral  ratios  over  the  low  frequency  (0-8  Hz)  and  high 
frequency  (6-12  Hz)  bands.  Note  the  reasonable  separation  for  the  low  frequency  band  and 
the  generally  poor  performance  over  the  high  frequency  band. 

We  applied  the  linear  discriminant  function  using  various  combinations  of  the  am¬ 
plitude  and  spectral  ratios  to  the  populations  of  8  earthquakes  and  8  mining  explosions. 
Table  1  shows  the  resiilts  using  the  holdout  classification  procedure.  That  is,  each  classifi¬ 
cation  is  made  using  a  learning  population  not  containing  the  event  to  be  classified.  The 
two  columns  are  the  number  of  misclassified  earthquakes  (as  explosions)  and  explosions  (as 
earthquakes)  respectively.  Dividing  the  entries  by  8  would  give  estimators  for  P(2|l)  and 
P(l|2)  respectively.  Results  using  the  minimum  discrimination  information  and  oi-entropy 
measures  described  in  the  next  section  are  also  included  for  purposes  of  comparison. 


Table  1:  Misclassifications  (Holdout) 


Method 

EQ 

EXP 

Amplitude  Ratio 

1 

1 

Amplitude-Low  Frequency  Spectral  Ratio 

1 

2 

Low  Frequency  and  High  Frequency  Spectral  Ratios 

0 

1 

Amplitude,  Low  and  High  Frequency  Spectral  Ratios 

0 

1 

Minimum  Discrimination  Information 

0 

1 

a- Entropy 

0 

0 

It  is  clear  that  the  amplitude  and  low  frequency  spectral  ratios  include  essentially  the 
same  information,  so  that  combining  them  does  not  work  well.  For  comparison  purposes, 
it  appears  that  the  log  amplitude  ratios  work  best  as  a  single  discriminator  and  that  the 
two  spectral  ratios  work  as  the  best  pair.  The  problem  event  is  always  the  first  explosion 
that  gets  classified  as  an  earthquake  in  the  holdout  procedure.  Using  amplitude  ratios,  the 
posterior  probability  that  the  first  explosion  is  an  earthquake  is  .88;  using  the  two  spectral 
ratios,  it  is  .97.  The  event  of  unknown  origin  at  Novaya  Zemyla  is  classified  as  an  explosion 
with  high  probability  (.98,  1.00)  by  the  two  procedures. 

1.4  Waveform  Discrimination 

For  the  waveform  approach  to  classification  we  regard  the  entire  time  realization  as  the 
signal  of  interest.  Because  of  the  separate  P  and  S  phases  that  always  seem  to  be  present,  it 
is  convenient  to  think  of  them  as  a  bivariate  process,  so  that  there  are  two  signal  vectors. 
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each  composed  of  a  window  (25  seconds  here)  of  one  of  the  phases.  It  is  fairly  evident 
that,  in  the  case  of  regional  events,  the  signals  may  be  regarded  as  zero  mean  processes 
and  that  the  differences  between  signals  from  earthquakes  and  explosions  might  be  due 
exclusively  to  differences  between  the  spectra.  The  two  P  and  S  signals  are  uncorrelated 
and  incoherent  (this  was  checked  for  the  sample  of  8  earthquakes  and  explosions  in  this 
study)  and  hence  can  be  regarded  as  imcorrelated  processes  with  unequal  P  and  S  spectra 
for  earthquakes  and  explosions. 

In  order  to  apply  the  quadratic  discriminant  ftmction,  which  would  be  appropriate 
for  distinguishing  between  the  earthquake  and  explosion  processes,  we  need  a  method  for 
computing  the  discriminant  function  (1.2).  Long  signals  make  this  computation  difficult 
and  it  is  conventional  to  apply  an  approximate  form  due  to  Whittle  [22],  say 

where  we  replace  •  by  P  or  S  depending  on  the  phase  to  be  considered,  X.(k)  is  the 
discrete  Fourier  transform  of  the  data  Xf.  and  fj-(i'k)  denotes  the  spectrum  for  phase  • 
under  hypothesis  Hj.  The  frequencies  are  of  the  form  Uk  =  k/T,  k  =  0, . . .  ,T  —  1.  Liggett 
(1971)  established  that  the  scaled  (by  1/T)  difference  between  (1.3)  and  (1.2)  converges 
almost  surely  to  zero.  The  optimal  statistic  for  testing  whether  the  sampled  bivariate 
series  is  from  Pi:  Earthquake  or  from  P2:Explosion  is  given  by 

Q  =  {diP  —  d2p)  +  (dis  —  dzs)-,  (1-4) 

where  we  accept  Hi  if  Q  >  0  and  H2  if  Q  <  0.  The  procedure  for  discrimination  using  a 
single  waveform  and  (3)  is  by  now  well  established  in  the  statistical  literature,  having  been 
applied  by  Alagon  (1989),  Cavanaugh  et  al  (1993),  Dargahi-Noubaryand  Laycock  (1981), 
Shumway  and  Unger  (1974),  Shumway  and  Blandford  (1974)  and  Shumway  (1982),  (1988) 
to  seismic  signal  discrimination  problems. 

It  is  useful  to  give  an  information  theoretic  interpretation  to  the  test  given  above  since 
this  exercise  leads  to  two  new  test  procedures  that  turn  out  to  be  robust  to  departures 
from  assiunptions.  KuUback  (1959)  has  developed  the  minimum  discrimination  information 
(MDI)  criterion  as  a  means  for  classifying  a  new  observation  into  Hi  or  H2-  Pinsker  (1964) 
derived  the  limiting  form  of  the  discrimination  information  under  the  assumption  that  one 
has  two  Gaussian  processes  differing  only  in  the  spectra.  For  the  MDI  criterion,  one 
compares  the  discrepancy  of  a  spectral  estimator  computed  from  the  sample  realization 
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Xt.,  say  /T-(i^fc),  with  fi.{uk)  and  f2.{i'k)  using 


I{fT;U)  =  I 


fj-M 


-log 


fr-iyk) 

fj-M 


(1.5) 


The  above  can  be  regarded  as  a  measure  of  the  match  between  the  sample  spectrum  and 
the  theoretical  value  fj-{i>k)  for  3  —  =  Pi^  earthquakes  and  j  =  2^  -  —  P,S  explosions. 

A  reasonable  procedure  would  be  to  look  at  the  difference 


=  i(fT  ,M  -  (1-6) 


between  the  matches  to  the  two  theoretical  spectra.  Since  we  want  the  discrepancy  between 
the  sample  spectrum  and  the  true  density  to  be  minimized,  it  is  clear  that  we  should  accept 
Hi  when  /(/i., /2-; /t-)  >  0  and  accept  H2  otherwise.  In  terms  of  the  overall  criterion, 
expressed  in  terms  of  both  phases  we  would  choose  Hi  (earthquake  )  when  the  sum  of 
(6)  over  the  ■  =  S,  •  =  P  exceeds  0  and  accept  H2  (explosion)  otherwise.  Note  that  for 
fT-{k>k)  =  |X(fc)p,  the  periodogram  estimator,  the  above  criterion  reduces  exactly  to  the 
quadratic  criterion  defined  in  (1.3).  Zhang  and  Taniguchi  (1993)  ,  (1994)  have  shown  the 
asymptotic  normality  of  the  MDI  criterion  and  that  the  misclassification  errors  converge 
to  zero.  They  have  also  shown  that  the  criterion  is  robust  to  departures  from  normality. 

Zhang  and  Taniguchi  have  also  suggested  the  Renyi  a-entropy,  (0  <  a  <  1)  see  Renyi 
(1961)  as  an  alternative  and  show  that  it  is  robust  to  both  non-Gaussian  departures  and 
peak  contamination.  Under  this  suggestion, define 

/r.)  =  2  E  {'°«(l  -  “  + 

and  accept  Hi  when 

^a(/l.,/2.;/T.)>0, 

where 

=  ea(/2.,/T-)  -  e^ifi-Jr-)  (1-8) 

In  terms  of  the  overall  criterion  involving  both  phases,  we  would  again  accept  Hi  when 
the  sum  of  (8)  over  ■  =  P,  -  =  S  exceeds  0. 

In  order  to  apply  the  discriminant  functions  defined  above,  we  need  to  have  estimators 
for  the  earthquake  and  explosion  spectra,  say  fi.(^)  /2-(z^).  These  can  be  taken  as 

predefined  values  if  no  training  sample  is  available  or  as  the  averages  of  the  earthquake  and 
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explosion  spectra  respectively  if  a  training  sample  is  available.  We  take  the  values  here  of 
the  average  earthquake  and  explosion  spectra  (see  Cavanaugh  et  al,  1993).  The  spectra 
were  computed  for  each  series  (no  taper)  over  a  fairly  broad  band  (2  Hz)  and  then  averaged 
separately  for  earthquakes  and  explosions.  The  event  to  be  classified  was  held  out  of  the 
averages.  Note  that  the  P  and  S  components  were  scaled  by  dividing  by  the  maximum  of 
the  P  component.  Figure  3  shows  the  spectra  average  earthquake  and  explosion  P  and  S 
spectra  for  a  fairly  broad  (2  Hz)  bandwidth  and  for  a  somewhat  more  narrow  bandwidth 
(.8  Hz).  Note  the  nulls  in  the  average  explosion  spectra  indicating  the  possibility  of  delay- 
fired  phenomena.  We  investigate  this  further  in  Part  2.  For  the  quadratic  and  information 
theoretic  detectors,  small  values  of  the  theoretical  spectra  can  cause  potential  distortions, 
so  several  cutoff  frequencies  (8  Hz,  12  Hz)  and  several  bandwidths  (.1  Hz,  .8  Hz,  2  Hz)  in 
(1.5)  and  (1.7)  were  tried;  overall  best  performance  seemed  to  be  attained  with  a  cutoff  of 
about  8  Hz.  The  discriminants  performed  equally  well  over  all  three  bandwidths  and  we 
chose  2  Hz  as  a  final  bandwidth. 

Table  1  shows  one  misclassification  (the  first  explosion  again)  for  the  discrimination 
information  and  no  misclassifications  for  the  a-entropy  with  a  =  .7.  Furthermore,  the 
a-entropy  shows  excellent  separation  between  the  earthquake  and  explosion  population 
as  shown  in  Figure  4.  The  minimum  discrimination  information  based  on  (1.5)  and  the 
optimal  quadratic  detector  based  on  (1.7)  both  separated  the  populations  as  well  but  there 
was  again  one  misclassification.  Overall,  the  a-entropy  yielded  the  largest  separation  of 
the  two  populations.  Again,  we  see  that  the  Novaya  Zemyla  event  falls  well  within  the 
explosion  group. 

1.5  Discussion 

In  Part  1,  we  have  reviewed  statistical  procedures  for  discriminating  between  earthquakes 
and  explosions  that  ensure  that,  for  a  set  explosion  false  alarm  rate,  the  detection  proba- 
bihty  will  be  maximized.  These  optimal  statistical  classification  procedures  turn  out  to  be 
linear  if  extracted  features  based  on  amplitude  and  spectral  characteristics  of  the  P  and  S 
phases  are  used.  Furthermore,  the  dimensionality  of  the  classification  function  is  usually 
relatively  low.  In  the  case  of  waveform  discrimination,  the  dimensionality  is  rather  high 
but  the  optimal  procedures  again  turn  out  to  be  fimctionally  dependent  on  the  spectra  of 
the  P  and  S  phases.  Optimal  discrimination  procedures  are  therefore  quadratic  functions 
of  the  original  waveforms  and  they  reduce  to  computing  functions  that  match  these  original 
waveforms  to  the  group  average  using  matching  functions  like  (1.3),  (1.5)  and  (1.7). 
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Figure  4:  Separation  Achieved  for  Information  Theoretic  Spectral  Matching  Measures.  The 

top  panel  shows  holdout  values  for  the  discrimination  information  whereas  the 
bottom  panel  shows  the  values  for  the  o'-entropy. 
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An  important  consideration  in  any  statistical  discrimination  procedure  is  the  extent  to 
which  the  particular  features  extracted  or  the  signal  waveforms  tend  to  conform  or  fail  to 
conform  to  underlying  assumptions  like  multivariate  normality  that  drive  the  theoretical 
derivations.  Hence,  for  extracted  features  from  the  particular  sample  of  Scandinavian 
earthquakes  and  explosions,  we  used  logarithms  and  took  ratios  of  components  of  the  two 
separate  P  and  S  waveforms.  For  the  waveform  discrimination  procedures,  theoretical 
results  of  Zhang  and  Taniguchi  (1993),  (1994),  available  for  the  discrimination  information 
and  for  the  a-entropy,  give  credence  to  the  claim  that  their  overall  performance  is  robust 
to  departures  from  normality  and  to  spectral  peak  contamination. 

In  the  particular  small  sample  of  8  earthquakes  and  8  explosions  and  the  unknown 
event  from  Novaya  Zemyla,  both  the  feature  extraction  and  the  waveform  discrimination 
approaches  performed  well  and  all  approaches  classified  the  event  from  Novaya  Zemyla 
with  the  mining  explosion  group.  All  methods  except  the  a-entropy  misclassified  the  first 
explosion  as  an  earthquake.  This  decision  was  a  strong  preference  of  the  featme  extraction 
method  but  was  on  the  borderline  for  the  quadratic  detector  and  for  the  discrimination 
information  measures.  The  a-entropy  still  classified  all  events  correctly  in  the  holdout 
sample.  These  latter  two  observations  tend  to  support  the  superiority  of  the  optimal 
statistical  classification  methods  as  applied  to  the  original  waveforms. 
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PART  2:  ANALYSIS  OF  RIPPLE  FIRED-SEISMIC  SIGNALS 


2.1  Echoes  and  Seismic  Signals 

Regional  seismic  monitoring  and  discrimination  capabilities  that  are  desirable  under 
a  potential  Comprehensive  Test  Ban  Treaty  (CTBT)  can  be  improved  by  developing  algo¬ 
rithms  and  new  procedures  for  distinguishing  between  earthquakes,  nuclear  explosions  and 
mining  explosions  of  various  kinds.  Much  effort  in  past  discrimination  studies  has  concen¬ 
trated  on  extracting  various  features  of  the  spectrum  that  are  characteristic  of  earthquakes, 
nuclear  explosions  or  mine  blasts. 

One  particular  spectral  feature  that  characterizes  some  mining  explosions  is  a  modu¬ 
lation  of  the  spectrum  introduced  by  a  ripple-fired  explosion.  A  ripple-fired  event  usually 
involves  detonation  of  a  number  of  explosions  that  are  often  regularly  grouped  in  space 
and  time.  Such  explosions,  known  as  quarry  blasts,  have  low  magnitudes  that  may  be 
close  to  those  of  nuclear  explosions  that  one  might  monitor  under  the  CTBT.  A  number  of 
authors  have  examined  various  aspects  of  this  problem  and  have  proposed  techniques  for 
analyzing  these  ripple-fired  seismic  signals.  Bamngardt  and  Ziegler  (1988)  have  used  aver¬ 
aged  spectra  and  cepstra  (see  Bogert  et  al,  1963)  to  look  for  ntdls  that  are  characteristic  of 
ripple-fired  events  at  regional  distances.  Hedlin  et  al  (1990)  noted  that  periodic  modula¬ 
tions  persisted  through  the  coda  of  ripple-fired  events  and  used  the  time  varying  spectrum 
or  spectogram  as  an  analysis  tool.  They  also  introduced  an  automated  procedure  based  on 
the  empirical  distribution  of  the  maximum  of  the  time  varying  cepstrum.  Chapman  et  al 
(1992)  look  at  spectral  modulation  in  the  time  varying  spectrum  or  sonogram  introduced 
by  some  specific  shot  geometries  and  deconvolve  using  the  complex  cepstrum.  Dysart  and 
PuUi  (1990)  use  the  cepstrum  as  a  measiure  of  complexity  for  discrimination,  noting  that 
explosions  generally  have  a  more  complicated  cepstrum,  whereas  earthquakes  have  a  less 
complex  form. 

The  above  approaches  rely  generally  on  spectral  nulls  for  qualitative  guidance  as  to 
whether  an  event  is  an  earthquake,  nuclear  explosion  or  mine  blast.  In  quarry  blasts, 
such  spectral  nulls  are  introduced  by  ripple-firing,  with  the  observed  series  composed  of 
a  signal,  added  and  delayed  by  the  sequence  of  shots,  with  an  additive  noise  component. 
In  this  paper,  we  will  take  a  new  look  at  the  predicted  patterns  in  frequency  and  time 
induced  by  the  signal  and  noise  model  for  ripple-fired  events.  In  the  frequency  domain, 
the  spectrum,  log  spectrum  and  cepstrum  are  well  known  measures  for  assessing  ripple 
firing.  The  spectrum  will  have  nulls  that  define  frequencies  corresponding  to  the  spacing 
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of  the  pulses  and  the  dtiration  of  the  pulse  train.  The  log  spectriun  enhances  the  resulting 
frequency  content  over  regions  of  low  power  and  the  cepstnnn  produces  peaks  at  the 
frequencies  present  in  the  log  spectrum,  enabling  one  to  estimate  spacing  delays  between 
successive  pulses  and  sometimes,  even  the  duration  of  the  sequence. 

A  new  time  domain  method  for  identifying  the  duration  spacing  of  the  ripple  fired 
pulses  is  introduced  which  uses  autoregressive  integrated  moving  average  (ARIMA)  mod¬ 
eling  in  conjimction  with  the  autocorrelation  and  partial  autocorrelation  functions  (ACF 
and  PACF)  to  estimate  the  delay  mechanism.  After  fitting  the  underlying  series  with  a 
low-order  autoregressive  model,  it  is  shown  that  the  ripple-fired  model  with  equal  delays 
will  induce  in  the  residuals  a  peak  in  the  PACF  at  the  spacing  of  the  ptdses  and  a  peak  in 
the  ACF  residuals  from  the  spacing  pulse  model  at  the  duration  value. 

The  frequency  domain  approach  again  uses  the  residuals  from  a  low-order  autoregres¬ 
sive  model  to  isolate  peaks  in  the  cepstrum  of  the  residuals.  The  main  contribution  here 
over  previous  methods  is  to  look  for  periodicity  in  the  log  spectrum  using  residuals  from 
previously  fitted  models  rather  than  the  underlying  series  where  such  peaks  are  generally 
obscured  by  the  frequency  response  of  the  earth-seismometer  system.  Stripping  the  sea¬ 
sonal  autoregressive  components  due  to  the  spacing  scalloping  uncovers  adjacent  spacing 
at  close  delays  that  were  obscured  in  the  original  spectrum. 

Both  the  time  domain  and  frequency  domain  methods  use  multiplicative  models  to  ap¬ 
proximate  the  spectrum  of  an  observed  process;  that  is,  we  assume  an  underlying  spectrum, 
fitted  by  a  autoregressive  series  of  low  order,  multiplied  by  seasonal  spectra  corresponding 
to  the  spacing  and  duration  parameters.  Both  methods  are  then  tested  on  two  simulated 
explosions  with  known  patterns  of  time  delays.  The  first  pattern  is  completely  regular  in 
spacing  whereas  the  second  pattern  has  some  variability  in  the  spacings.  Both  methods 
work  well  in  this  case.  Applying  the  methods  to  two  Scandinavian  explosions  and  two 
earthquakes  leads  to  predicted  spacings  for  the  the  explosions  and  to  nonsignificant  peaks 
in  the  earthquake  population. 

2.2  Frequency  and  Time  Domain  Methods 

The  general  model  for  an  observed  series  j/t  that  contains  echoes  can  be  written  in 
terms  of  a  train  of  observed  signals  that  have  been  delayed  and  added,  plus  a  noise,  i.e., 

n— 1 

j/t  =  ^  ajSt-Ti  +  (2-1) 

j=o 
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where  rit  is  additive  noise,  To,  Ti , . . . ,  Tn-i  are  the  time  delays  associated  with  the  n  arrivals 
generated  by  the  ripple  firing  and  oq,  <^2,  •  -  • ,  <3n-i  are  the  amplitudes  of  the  arrivals.  Often, 
by  convention,  oq  =  1  and  To  =  0  so  that  the  initial  signal  appears  with  no  delay  and  unit 
amplitude.  A  similar  problem  that  has  occurred  in  the  past  is  testing  for  P-pP  delays  in 
long  period  data  as  by  Shiunway  and  Blandford  (1978)  who  used  a  version  of  (2.1)  with 
n  =  2,  oo  =  IjTo  =  0  to  estimate  the  parameters  Oi,  the  amplitude  of  the  reflection  and 
Ti,  the  P-pP  delay  time.  The  predicted  spectrum  Paii^)  was  derived  from  the  Von  Seggem 
and  Blandford  (1972)  source  theory.  Earlier  work  of  Hannan  and  Thomson  (1974)  used  a 
likelihood  approach  with  frequency  dependent  delay  times. 

Many  kinds  of  multipath  reflections,  etc.,  can  produce  the  behavior  predicted  by  the 
signal  model  (2.1)  and  we  shall  be  interested  in  the  patterns  that  might  be  induced  by 
regular  ripple-firing  of  the  kind  expected  with  mining  explosions.  As  an  example,  taking 
Tj  =  jd  and  aj  =  1  in  the  model  implies  that  we  are  interested  primarily  in  a  sequence 
of  pulses  of  equal  amplitude,  occurring  at  0,  d,  2d, . . . ,  (n  —  l)d.  This  reduces  (2.1)  to  the 
form 

n— 1 

>=0 

where  the  delays  are  at  multiples  of  an  underlying  delay  d  and  the  amplitudes  are  the 
same.  This  implies  that  there  are  n  pulses  spaced  at  interval  d  and  we  will  refer  to  d  as 
the  spacing  parameter  and  to  the  quantity  nd  as  the  duration  of  the  signal.  Hopefully, 
this  might  produce  a  good  correspondence  with  what  could  be  expected  from  a  ripple-fired 
pulse  train. 

In  this  case,  there  are  the  two  parameters  d  and  n  that  must  be  identified,  along  with 
the  spectra  of  the  signal  and  noise,  say  Ps(u})  and  Pn{i^),  where  u?  is  the  frequency  in 
cycles  per  unit  time.  Assuming  that  the  signal  and  noise  processes  are  independent,  the 
spectrum  of  the  output  yt  should  be  of  the  form 

P,{w)  =  \A(^)fP.(u)  +  P„(u,),  (2.3) 


where 


iA(c.)r  = 


2  sin^{Trujnd) 


(2.4) 


sin'^{wujd) 

The  above  shows  that  one  should  expect  the  underlying  signal  spectrum  to  be  multiplied 
by  a  frequency  response  function  that  is  periodic  with  frequencies  proportional  to  d  and 
nd.  Hence,  the  cepstrum,  which  is  the  spectrum  of  the  log  spectrum,  should  show  at  least 
two  peaks  at  delays  of  d  and  nd,  corresponding  to  the  spacing  and  duration  of  the  pulse 
train. 
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To  show  an  example,  consider  the  contrived  event  in  Figure  5  which  was  constructed 
using  a  simulated  2nd  order  autoregressive  signal  process,  where  a  order  process,  de¬ 
noted  by  AR(p)  in  this  paper,  is  of  the  form 

p 

St  =  <f>jSt-j  +  wt,  (2-5) 

i=i 

with  Wt  is  a  white  noise  input  process  with  variance  .  In  Figme  5,  we  have  taken 
p  =  2,  (f>i  =  l,(f>2  =  —  .6,  =  1,  and  the  process  has  a  smooth  spectrum  with  a  single  peak 

that  emvdates  that  of  an  explosion.  The  process  is  modulated  with  a  decaying  exponential 
function  of  the  form 

gt  =  (2.6), 

where  log^i  =  —7,^2  =  -01,  producing  the  signal  St  =  gtSt-  The  resulting  signal  is  then 
delayed  according  to  the  model  from  (2.2)  with  no  noise  and  n  =  5,  and  d  =  8.  Hence 
there  are  five  piilses,  each  separated  by  eight  time  points,  leading  to  a  signal  with  spacing 
5  points  and  duration  40  (5  X  8)  points.  The  time  response  of  the  pulses  is  shown  in  the 
upper  left  hand  panel  of  Figure  6.  The  waveform  on  the  lefthand  side  of  Figure  5  has  a 
form  that  emulates  many  explosions  such  as  those  found,  for  example,  in  Cavanaugh  et 
al  (1993).  In  the  right  hand  half  of  Figure  5  is  shown  a  process  that  has  been  delayed 
with  irregular  spacings;  for  example,  we  took  To  =  Ti  =  8,  T2  =  T3  —  10,  r4  =  8,  Ts  =  10. 
The  simulated  explosions  in  Figure  5  are  quite  similar  to  the  real  Scandinavian  explosions 
shown  in  Figure  11  and  we  assume  that  the  generating  procedure  is  reasonable.  We  assiime 
implicitly  a  sampling  rate  of  40  points  per  second,  leading  to  a  folding  frequency  of  20  Hz. 

Bogert  et  al  (1962)  proposed  computing  the  cepstrum,  i.e.  the  spectrum  of  the  log 
spectrum,  and  looking  for  peaks  that  corresponded  to  periodicities  in  the  log  spectrum. 
The  peaks  might  be  identified  with  the  spacing  and  duration  of  the  sequence  of  pulses, 
since  the  logarithm  of  the  product  of  the  signal  spectrum  and  the  modulating  function  in 
(2.3)  will  be  roughly  proportional  to 

log  +  2  log  sin{'KU)nd)  —  2  log  sin{'Kud)  (2-7) 

which  displays  it  as  the  sum  of  the  signal  spectrum  and  two  nonlinear  periodic  functions. 
To  get  an  idea  of  what  this  should  look  like,  we  have  computed  the  cepstrum  of  the  regular 
pulse  train  in  Figure  6  and  the  result  shows  major  peaks  at  delays  of  8  and  40  points  with 
intermediate  values  at  multiples  of  the  spacing  interval  of  8  points.  Hence  the  spacing  (8 
points)  and  duration  (8  x  5  =  40  points)  appear  as  the  major  components  of  the  cepstrum. 
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Figure  5:  Simulated  Series  and  2nd-0rder  AR(2)  Autoregressive  Spectra.  The  left  panel 
shows  a  simulated  series  formed  from  an  echo  model  with  n  =  5  replicates  of  an 
AR(2)  series  spaced  at  fixed  delays  with  spacing  d  =  8.  The  right  panel  shows 
the  same  series  with  n  =  6  pulses,  2  at  d  =  8  followed  by  2  at  d  =  10  followed  by 
1  at  d  =  8  and  1  at  delay  d  =  10.  Spectra  measured  in  cycles  per  point  (.5  cycles 
per  point=20  Hz). 
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PACF  of  Pulse  Train  (d=8,  n= 


Figure  6:  Theoretical  Pulse  Train  Corresponding  to  the  Simulated  Series  in  the  Left  Panel 
of  Figure  2.1.  Also  shown  are  the  cepstrum,  the  PACF  of  the  pulse  train  and 
the  ACF  of  the  seasonally  differenced  pulse  train.  Note  the  well  defined  peaks 
corresponding  to  d  =  8  and  n  =  5  (spacing  and  duration  of  8  and  40  points). 


Significance  values  for  peaks  in  the  cepstnim  can  be  approximated  using  the  fact  that 
the  unsmoothed  cepstral  values  are  approximately  distributed  as  a  chi-squared  random 
variable  with  2  degrees  of  freedom.  This  is  argued  from  the  fact  that  logarithms  of  adja¬ 
cent  discrete  raw  spectra  can  be  regarded  as  approximately  Gaussian  and  independent  so 
that  the  Fourier  transforms  of  the  same  axe  also  Gaussian.  Hence,  the  cepstrum,  as  the 
squared  Fourier  transform,  should  have  approximately  a  chi-squared  distribution.  Then, 
the  behavior  of  the  peak  values  can  be  compared  using  the  fact  that  a  chi-squared  random 
variable  with  2  degrees  of  freedom  would  be  expected  to  exceed  3  times  its  theoretical 
value  only  about  5%  of  the  time.  Hence,  peahs  that  emerge  from  the  underlying  cepstral 
pattern  can  be  assumed  to  be  significant  if  they  are  more  than  three  times  the  underlying 
values. 

An  alternative  time  domain  method  that  does  not  appear  to  have  been  tried  is  based  on 
converting  the  signal  model  to  an  approximately  equivalent  time  domain  ARIMA  process 
and  then  applying  the  usual  diagnostic  techniques.  Writing  the  simplified  model  (2.2)  in 
terms  of  the  usual  shift  operator  Bst  =  Sf-i  leads  to 


n— 1 

yt  =  ^  +  ut 

i=0 

-  (1  _  (2-8) 

which  exhibits  the  signal  part  of  the  model  as  a  multiplicative  model  with  a  seasonal 
autoregressive  part  with  period  d  and  a  seasonal  moving  average  part  with  period  nd.  It 
follows  that  we  might  look  for  a  ratio  of  seasonal  polynomials  of  the  form 


a 


{B)  = 


(1  - 


where  $  and  0  vary  but  may  be  close  to  one.  If  we  assume  that  the 
can  be  satisfactorily  fitted  with  a  simple  autoregressive  model  of  order 
express  the  overall  model  in  the  multiplicative  form 


(2.9) 

imderlying  signal 
p,  then  we  might 


where 


yt 


(B) 

<t>{B) 


<f>{B)  =  1-^2 
j-i 


(2.10) 


(2.11) 
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is  the  autoregressive  model  for  the  modulated  signal  and  Wt  is  a  white  noise  process  with 
variance  In  general,  we  ignore  the  nonlinear  modulating  function  gt  and  concentrate  on 
fitting  a  low  order  autoregressive  model  to  the  imderlying  signal  first  so  that  we  estimate 
the  polynomial  <f>{B)  We  take  the  order  to  be  that  suggested  by  the  partial  autocorrelation 
function  (PACF)  which  tends  to  be  zero  after  the  given  order.  Note  that  one  should  be 
careful  not  to  fit  too  high  an  order  so  as  to  obscure  components  of  the  polynomial  cx(B). 

Now  if  one  follows  the  usual  diagnostic  procedure  (see,  for  example,  Shumway,  1988), 
we  would  look  at  residuals  from  the  fitted  AR(p),  where  we  should  see  a  pattern  suggesting 
the  form  of  Oi(B).  The  following  pattern  should  be  present.  The  PACF  should  have  a  peak 
at  d  whereas  the  ACF  should  have  a  peak  at  nd.  As  a  practical  matter,  it  works  better  to 
look  for  the  peak  in  the  PACF  at  d  and  then  to  fit  some  seasonal  AR  model;  the  ACF  of 
the  residuals  from  this  model,  say  (1  —  ^B^)yt  tend  to  show  peaks  at  nd.  Applying  this 
procedure  to  the  pulse  train  in  Figme  6,  for  example,  shows  a  well-defined  peak  at  d  =  8 
points  in  the  PACF  of  original  series  and  a  well-defined  peak  at  nd  =  40  points  in  the  ACF 
of  the  residuals  (1  —  B^)yt.  It  follows  that  the  PACF  and  ACF  of  the  seasonal  residuals 
match  the  pattern  that  would  be  expected  for  a  pulse  train.  The  cepstrum  of  the  series 
of  pulses  shows  primary  peaks  at  d  =  8  points  and  nd  =  40  points  as  would  be  expected 
from  the  underlying  model  which  is  5  pulses  spaced  by  8  points.  Note  also  the  peaks  in 
the  cepstrum  at  lags  8, 16  and  32  points. 

We  check  the  diagnostic  capabilities  of  the  cepstrum  and  the  time  domain  method 
by  applying  the  same  technique  to  the  simulated  explosions  in  Figure  5.  In  this  case,  the 
following  diagnostic  procedure  was  followed  for  the  first  simulated  signal. 

(i)  Use  the  fact  that  the  PACF  in  Figure  7  shows  peaks  at  the  first  two  lags  and  then  a 
break,  indicating  that  the  underlying  spectrum  of  the  process  might  be  satisfactorily 
emulated  by  fitting  the  second-order  autoregressive  model  of  the  form 

(1  -  <i>2B  -  =  Wt  (2.12) 

(ii)  The  partial  autocorrelation  function  (PACF)  of  the  estimated  residuals  Wf  =  ^(B)yt 
from  the  above  model  shows  a  strong  peak  at  lag  8,  indicating  the  denominator  of 
(2.9)  holds  with  d  =  8  The  model  now  includes  a  seasonal  autoregressive  component, 
leading  to 

(1  -  $iR»)(l  -  <f>2B  -  <i>2B^)yt  =  Wf  (2.13) 

as  the  current  proposed  model.  Note  that  we  might  also  have  included  a  seasoned 
difference  term  with  =  1  but  some  experimentation  shows  that  the  estimated 
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Figure  7:  Autocorrelation  Functions  (ACF)  and  Partial  Autocorrelation  Functions  (PACF) 
of  the  First  Simulated  Series.  Top  panel  shows  ACF  and  PACF  of  the  original 
series;  the  middle  panel  shows  ACF  and  PACF  from  residuals  of  the  AR{2)  model; 
the  bottom  panel  shows  ACF  and  PACF  from  residuals  of  SAR(8)xAR(2)  model. 
Lags  are  in  points  (1  point=.025  sec.). 


values  for  contrived  data  are  not  that  near  the  unit  circle.  Furthermore  the  ACF  of  the 
residuals  does  not  exhibit  the  required  peaks  at  d,  2d,  3d, . . .  that  would  be  expected 
for  the  nonstationary  seasonal  model. 

(iii)  Autocorrelations  (ACF)  and  partial  autocorrelations  (PACF)  for  residuals  from  the 
model  above  are  shown  in  the  bottom  panel  of  Figure  7  and  they  have  the  required  peak 
at  lag  40,  indicating  that  n  =  5,  d  =  8  are  the  required  values  for  the  parameters.  The 
final  model  then  has  the  seasonal  moving  average  component  with  nd  =  40,  leading 
to  the  model 

(1  -  -  hB  -  <t>2B^)yt  =  (1  -  (2.14) 

This  is  not  exactly  the  Box-Jenkins  (see  Shumway,  1988)  form  because  the  seasonalities 
on  the  right  hand  and  left  hand  side  do  not  match. 

We  note  that  the  results  for  the  cepstrum  of  the  residuals  Wt  =  (f)(B)yt  ,  shown  in 
Figure  8,  also  match  almost  exactly  the  expected  forms  displayed  in  Figure  6.  The  only 
difference  in  the  time  domain  results  is  some  additional  noise  in  the  ACF  and  PACF  but 
the  pattern  that  leads  to  concluding  that  the  polynomial  in  (2.8)  is  present  still  persists. 
We  conclude  that  the  pattern  generated  by  a  fixed  sequence  of  shots,  detonated  at  equal 
spacings  can  be  detected  by  either  the  time  or  frequency  domain  methods  in  this  case. 

For  the  second  simtilated  explosion,  containing  irregular  delays  at  ten  points  (3  de¬ 
lays)  and  at  eight  points  (3  delays),  the  situation,  shown  in  Figure  9,  is  somewhat  more 
complicated.  We  still  see  the  peaks  at  lags  1  and  2  in  the  PACF,  suggesting  fitting  first 
an  AR(2)  model  as  before.  Residuals  from  this  model,  shown  in  the  middle  panel  of  Fig¬ 
ure  9,  suggest  that  there  are  delays  present  at  8  and  10  points,  using  the  PACF.  Now, 
fitting  a  seasonal  5Ai2(10)  in  the  strongest  peak  leads  to  residuals  supporting  the  peak  at 
8  points,  indicating  a  second  kind  of  delay.  Residuals  from  the  model  with  both  seasonal 
peaks  (1  —  $2-B^)(1  —  ^iB^°)yt  do  not  show  weak  peaks  at  any  particular  duration,  say 
3x10 -1-3x8  =  54  as  would  be  expected.  The  cepstriun,  shown  in  Figure  10,  shows  double 
peaks  at  8  and  10  points  as  well  as  a  peak  at  lag  54,  corresponding  to  the  duration.  There 
is  an  additional  peak  at  the  intermediate  value  corresponding  to  28  points.  Computing 
residuals  of  the  form  (1  —  ^iB^)(l  —  ^iB  —  ^2B^)yt  still  shows  the  peak  at  10,  28  and  55 
points. 

In  the  next  section,  we  apply  the  two  single-channel  single-event  techniques  above  to 
two  earthquakes  and  explosions  from  the  Scandinavian  regional  event  population  consid¬ 
ered  by  Blandford  (1993)  and  by  Cavanaugh  et  al  (1994). 
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delay 


Figure  8:  Cepstral  Analysis  of  the  First  Simulated  Series.  Top  panel  shows  cepstrum  of  the 
AR(2)  residual  series;  bottom  panel  shows  cepstrum  of  SAR(S)  x  AR(2)  residual 
series.  Lags  are  in  points  (1  point=.025  sec.). 
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Figure  9;  Autocorrelation  Functions  (ACF)  and  Partial  Autocorrelation  Functions  (PACF) 
of  the  Second  Simulated  Series.  Top  panel  shows  ACF  and  PACF  of  the  original 
series;  the  middle  panel  shows  ACF  and  PACF  from  residuals  of  the  AR{2)  model; 
the  bottom  panel  shows  ACF  and  PACF  from  residuals  of  5'A.i?(10)  x  AR{2) 
model.  Lags  are  in  points  (1  point=.025  sec.). 
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Figure  10:  Cepstral  Analysis  of  the  Second  Simulated  Series.  Top  panel  shows  cepstrum 
of  the  AR(2)  residual  series;  bottom  panel  shows  cepstrum  of  SAR{8)  x  AR{2) 
residual  series.  Lags  are  in  points  (1  point=.025  sec.). 
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2.3  Analysis  of  the  Scandinavian  Events 

In  order  to  get  a  sense  for  the  structures  of  several  typical  earthquakes  and  explosions 
in  the  Scandinavian  popiilation,  consider  the  two  earthquakes  and  two  explosions  taken 
from  those  considered  by  Blandford  (1993).  The  events  chosen  were  two  explosions  (local 
magnitudes  2.32  and  2.59,  taken  from  Table  1  of  Blandford)  and  two  earthquakes  (local 
magnitudes  3.26  and  4.40  taken  from  Table  2  of  Blandford).  All  events  chosen  by  Blandford 
were  on  or  near  land  and  were  distributed  uniformly  over  Scandinavia  to  minimize  the 
possibility  that  discriminators  might  be  keying  on  location  or  land-sea  dilferences.  The 
P-phases  from  the  two  explosions  and  the  two  earthquakes  are  plotted  in  Figures  11  and 
16  respectively.  Note  that  the  waveform  of  the  second  explosion,  shown  in  the  right  panel 
of  Figure  11,  resembles  the  contrived  events  shown  in  Figure  5.  We  consider  analyzing  the 
two  explosions,  labeled  5P  and  6P  and  the  two  earthquakes,  labeled  5P  and  6P  by  the 
methods  proposed  in  the  previous  section. 

(i)  For  the  first  explosion,  labeled  EXP  5P  in  Figxure  11,  it  is  necessary  to  first  accoimt  for 
the  primary  contributor  to  the  smooth  peak  part  of  the  spectrum.  Figure  12  shows  the 
ACF  and  PACF  of  the  original  P-phase  and  it  is  clear  from  the  PACF  that  the  order 
of  the  vmderlying  process  could  be  taken  as  p  =  4.  Fitting  an  autoregressive  AR{A) 
model  leads  to  the  peak  spectrum  shown  in  Figtire  11.  The  residuals  from  this  model, 
given  in  Figure  12,  have  ACF  and  PACF  showing  peaks  at  delays  d  —  17, 19  points 
indicating  reflections  occmrring  at  these  two  spacings.  Fitting  a  seasonal  5Ai2(19) 
model  to  the  AR{A)  residuals  shows  that  the  peak  at  17  persists  and  we  conclude 
tentatively  that  there  is  activity  at  spacings  of  17  and  19  points,  corresponding  to 
about  .5  seconds.  The  cepstral  analysis  in  Figure  13  yields  a  similar  conclusion. 
There  axe  peaks  at  lag  19  and  132  points  in  the  residual  AR{A)  cepstrum,  suggesting 
that  a  seasonal  Ai2(19)  model  be  tried  and  that  the  dmration  of  the  pulse  sequence 
is  about  132  points  or  3.3  seconds  with  a  spacing  of  about  .5  seconds.  The  duration 
should  also  appear  in  the  PACF  at  lag  132  but  we  did  not  compute  this  due  to  the 
length  of  lag  required.  Fitting  the  seasonal  SAR{19)  model  xmcovers  a  second  peak 
at  17  points  in  the  AR{A)  x  5A72(19)  residuals,  suggesting  that  there  are  two  spacings 
operating  in  the  sequence  as  with  the  simulated  data  analyzed  in  Figure  10.  We 
conclude  tentatively  that  there  might  be  a  pulse  train  of  duration  3.3  seconds  with 
pulses  spaced  irregularly  at  about  .5  seconds.  Note  that,  in  this  example,  the  ACF 
and  PACF  resolved  the  mixed  delays  on  the  initial  examination  of  the  AR(4)  residuals 
whereas  the  cepstral  analysis  needed  to  have  the  peak  at  lag  19  stripped  away. 
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Figure  11:  P-Pha^es  From  Two  Scandinavian  Mining  Explosions  and  ^  Orto  AH(4) 

Autoregressive  Spectra.  Sampling  rate  is  40  points  per  second.  Spectra  measured 

in  cycles  per  point  (.5  cycles  per  point=20  Hz). 
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Figure  12;  Autocorrelation  Functions  (ACF)  and  Partial  Autocorrelation  Functions  (PACF) 
of  the  First  Mining  Explosion.  Top  panel  shows  ACF  and  PACF  of  the  original 
series;  the  middle  panel  shows  ACF  and  PACF  from  residuals  of  the  AR(4)  model; 
the  bottom  panel  shows  ACF  and  PACF  from  residuals  of  5Ai2(19)  x  AR{2) 
model.  Lags  are  in  points  (1  point=.025  sec.) 
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Figure  13; 


Cepstral  Analysis  of  the  First  Mining  Explosion  .  Top  panel  shows  cepstrum  of 
the  AR(4)  residual  series;  bottom  panel  shows  cepstrum  of  SAR{19)  x  AR{2) 
residual  series.  Lags  are  in  points  (1  point=.025  sec.). 
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(ii)  The  second  explosion,  labeled  EXP  6P  in  Figure  11,  also  was  fitted  reasonably  well  by 
an  AR{4)  model  since  the  PACF  in  Figure  14  is  essentially  zero  after  lag  4.  Looking 
at  the  PACF  (and  ACF)  of  the  Ai2(4)  residuals  picks  up  a  peak  at  lag  12  and  we 
tentatively  identify  an  AR{4)  x  SAR{12)  model  for  the  next  step.  Residuals  from 
the  above  model  tend  to  show  a  peak  at  lag  24  in  the  PACF,  implying  a  SAR(24:) 
component  identified  with  the  duration  of  two  pulses,  spaced  at  lag  12  (.3  seconds). 
Looking  at  the  cepstral  analysis  in  Figure  15,  we  see  first  the  peak  at  12  in  the  Ai2(4) 
residuals  and  then  peaks  at  12,  24,  36,  42,  57  points  in  the  AR{A)  x  SAR{12)  residuals. 
This  suggests  the  possibility  of  5  pulses  spaced  at  about  .3  seconds  rather  than  the 
two  pulses  implied  by  the  analysis  of  the  ACF  and  PACF.  Here,  we  note  that  the 
cepstrum  seems  to  focus  better  on  the  pattern  for  the  higher  reflections. 

(iii)  The  first  earthquake,  EQ  5P,  shown  in  Figure  16,  has  substantial  coefficients  in  the 
PACF,  as  shown  in  Figure  17,  for  lags  up  to  4.  The  Ai?(4)  spectrum  shown  in  Figure 
16  has  a  rather  strong  low  frequency  ripple  around  .125  cycles  per  point  or  5  Hertz. 
Fitting  a  lower  order  AR{2)  to  the  underlying  spectrum  omitted  the  dip  but  did  not 
change  the  conclusions  below.  Figure  17  shows  the  ACF,  PACF  and  cepstrum  of 
the  AR{A)  residuals  and  we  do  not  see  any  singular  peaks  in  the  two  time  domain 
measures  or  any  prominent  peaks  in  the  cepstrum.  The  largest  peak  in  the  cepstrum 
at  delay  7  is  only  about  2  times  the  lowest  possible  baseline  level. 

(iv)  The  second  earthquake,  EQ  6P,  shown  in  the  right  hand  panel  of  Figure  16,  has  a 
PACF,  shown  in  Figure  18,  that  disappears  after  lag  6,  implying  that  an  AR{6)  model 
should  be  fitted  to  the  underlying  spectrum.  The  ACF,  PACF  and  cepstrum  of  the 
AR(6)  residuals  again  do  not  have  any  particularly  prominent  peaks,  indicating  again 
an  absence  of  reflections. 

To  conclude,  we  note  that  the  explosions  both  displayed  strong  indications  of  reflec¬ 
tions  at  equal  spacings  (Explosion  6P)  or  irregular  spacings  (Explosion  5P).  In  contrast, 
the  two  earthquakes  displayed  a  more  complex  structure  for  the  underlying  signal,  forcing 
us  to  fit  a  higher  order  autoregressive  model  for  that  component.  The  residuals,  however, 
had  a  more  simple  pattern  with  no  indications  of  reflections  that  might  be  characteristic 
of  a  ripple-fired  event. 

2.4  Discussion 

We  have  developed  a  frequency  domain  method  using  the  peaks  in  the  cepstrum  to 
identify  the  spacing  and  duration  of  echoes  in  a  simple  ripple-fired  model  and  shown  that 
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Figure  14:  Autocorrelation  Functions  (ACF)  and  Partial  Autocorrelation  Functions  (PACF) 
of  the  Second  Mining  Explosion.  The  top  panel  shows  ACF  and  PACF  of 
the  original  series;  the  middle  panel  shows  ACF  and  PACF  from  residuals  of 
the  AR{A)  model;  the  bottom  panel  shows  ACF  and  PACF  from  residuals  of 
SAR(19)  X  AiJ(4)  model.  Lags  are  in  points  (1  point=.025  sec.). 
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Figure  15: 


Cepstral  Analysis  of  the  Second  Mining  Explosion  .  The  top  panel 

of  the  AJJ(4)  residual  series;  the  bottom  panel  shows  cepstrum  of  X 

AR{A)  residual  series.  Lags  are  in  points  (1  point-.025  sec.). 
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Figure  16;  P-Phases  From  Two  Scandinavian  Earthquakes  and  Fourth  Order  ylii(4)  and 
Sixth  Order  AR{6)  Autoregressive  Spectra.  Sampling  rate  is  40  points  per  second. 
Spectra  measured  in  cycles  per  point  (.5  cycles  per  point=20  Hz). 
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Figure  17:  ACF,  PACF  and  Cepstra  of  the  First  Earthquake.  The  top  panel  shows  ACF 
and  PACF  of  the  original  series;  the  middle  panel  shows  ACF  and  PACF  from 
residuals  of  the  AR(4)  model;  the  bottom  panel  shows  the  cepstrum  computed 
from  residuals  of  AR{4)  model.  Lags  are  in  points  (1  point=.025  sec.). 
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Figure  18:  ACF ,  PACF  and  Cepstra  of  the  Second  Earthquake.  The  top  panel  shows  ACF 
and  P ACF  of  the  original  series;  the  middle  panel  shows  ACF  and  PACF  from 
residuals  of  the  Ai?(6)  model;  the  bottom  panel  shows  the  cepstrum  computed 
from  residuals  of  AR(6)  model.  Lags  are  in  points  (1  point=.025  sec.). 


39 


it  performs  satisfactorily  on  contrived  and  real  data.  A  new  time  domain  method  using  the 
ACF  and  PACF  to  build  a  seasonal  autoregressive  moving  average  model  for  the  ripple¬ 
firing  phenomenon  also  worked  well  on  the  simulated  and  real  data.  It  should  be  noted  that 
the  next  step  in  the  procedure  should  be  to  fit  more  closely  some  of  the  models  implied  by 
the  identification  procedures  developed  in  this  paper.  That  is,  the  diagnostic  procedures 
can  be  used  to  suggest  models  of  the  form  given  in  (2.2)-(2.4)  in  the  simplest  case.  We 
need  to  develop  methods  for  comparing  models  and  for  estimating  the  parameters  n  and 
d  by  maximum  likelihood  or  some  other  efficient  method. 


To  illustrate  one  such  possible  procedure,  let  Y'r{ujk)  denote  the  discrete  Fourier  trans¬ 
form  of  yti  evaluated  at  frequencies  of  the  form  =  27rk/T,  k  =  0, ...  ,T  —  1  and  write 
the  approximate  log  likelihood  function 


T-l  T-l 

log  i:(n,  d)  oc  -  ^  logPj,(a;ifc)  - 

Jfc=0  k=0 


\YtM\'‘ 
P,M  ■ 


(2,15) 


Maximizing  the  above  function  with  respect  to  n  and  d,  where  Py(uk)  is  given  by  (2.3) 
and  (2.4),  would  lead  to  the  best  fitting  simple  model.  The  approach  could  be  extended 
to  more  irregular  patterns  such  as  that  given  by  (2.1)  since  there  is  no  particular  reason 
why  the  amplitudes  and  time  delays  in  (2.1)  should  have  unit  amplitudes  and  be  equally 
spaced  as  is  implied  by  Equation  (2.2)  and  the  arguments  following  it.  Allowing  the 
amplitudes  aj,j  =  0, 1, . . . ,  n  —  1  and  time  delays  To, Ti, . . . , Tn-i  to  take  general  values 
involves  adding  additional  parameters  to  the  model.  If  general  amplitudes  and  time  delays 
are  allowed,  there  is  the  problem  of  distinguishing  between  ripple-fired  scalloping  and 
scalloping  introduced  by  multipath  effects  or  by  local  site  anomalies. 


The  above  complications  suggest  that  an  array  deconvolution  approach  also  might  be 
helpful.  Suppose  that  we  collect  a  number  of  explosion  signals,  all  recorded  at  the  same 
array,  so  that  we  observe  the  jth  signal  at  the  ith  array  as 


Vijt  —  fit  *  d"  f^ijt  (2.16) 

for  i  =  1, 2, . . . ,  r  arrays  observing  j  =  1, 2, . . . ,  m  signals.  The  notation  *  in  (2.16)  denotes 
convolution.  The  model  separates  effects  characteristic  of  the  site,  say  Vit,  the  signals 
Sjt.  Shumway  and  Der  (1984)  have  developed  a  method  for  simultaneously  estimating  the 
receiver  functions  by  maximmn  likelihood  and  then  deconvolving  the  signals.  As  a  part  of 
the  procedure,  the  imderlying  signal  spectra  Fs(w)  and  noise  spectra  P„(a;)  can  be  easily 
estimated  as  in  Der  et  al  (1992)  so  that  the  deconvolved  signals  will  presumably  contain 
only  the  common  reflective  pattern.  We  have  not  yet  investigated  any  of  the  powerful 
capabilities  of  this  array-based  methodology. 
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